#ifndef PCPS_OPERATOR_HEADER
#define PCPS_OPERATOR_HEADER

#include <src/pcps.h>

namespace pcps {

  //--------------------------------------------------------------------------------
  // pcps::OpBase -- base class for second quantized operators
  //--------------------------------------------------------------------------------

  template <class SCALAR> class OpBase {

    // constructors
    public:
      OpBase() {}
      OpBase(const int ns, const SCALAR coeff) : _sites(ns, -1), _offsets(ns, 0), _coeff(coeff) {}
      OpBase(const std::vector<int> & sites, const std::vector<int> & offsets, const SCALAR coeff)
      : _sites(sites), _offsets(offsets), _coeff(coeff) {}
      OpBase(const OpBase &other) : _sites(other._sites), _offsets(other._offsets), _coeff(other._coeff) {}

    // destructor (virtual so that this class's children's destructors will be called when base pointers are deleted)
    virtual ~OpBase() {}

    // assignment operator (deep copy)
    public:
      OpBase &operator=(const OpBase &other) {
        if (this == &other) return *this;
        _sites.assign(other._sites.begin(), other._sites.end());
        _offsets.assign(other._offsets.begin(), other._offsets.end());
        _coeff = other._coeff;
        return *this;
      }

    // smart pointers
    public:
      typedef boost::shared_ptr<OpBase> ptr;
      typedef boost::shared_ptr<const OpBase> const_ptr;

    // data access functions
    public:
      const std::vector<int> & sites() const { return _sites; }
      const std::vector<int> & offsets() const { return _offsets; }
      const SCALAR & coeff() const { return _coeff; }
      std::vector<int> & sites() { return _sites; }
      std::vector<int> & offsets() { return _offsets; }
      SCALAR & coeff() { return _coeff; }

    // other public functions
    public:
      virtual void get_rhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const = 0;
      virtual void get_uhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const = 0;
      virtual double action_value(const bool particles_are_spins,
                                  const bool particles_are_fermions,
                                  const int configs_per_site,
                                  const int * occs) const = 0;

    // data
    protected:
      std::vector<int> _sites;
      std::vector<int> _offsets;
      SCALAR _coeff;

  };

  //--------------------------------------------------------------------------------
  // pcps::OpNN -- class for two number operators acting on different sites
  //--------------------------------------------------------------------------------

  template<class SCALAR> class OpNN : public OpBase<SCALAR> {

    // constructors
    public:
      OpNN() : OpBase<SCALAR>(2, 0.0) {}
      OpNN(const int site0, const int site1, const SCALAR coeff)
        : OpBase<SCALAR>(2, coeff)
      {
        if (site0 == site1)
          throw pcps::Exception("The two sites of OpNN must be distinct");
        this->_sites.at(0) = site0;
        this->_sites.at(1) = site1;
        this->_offsets.at(0) = 0;
        this->_offsets.at(1) = 0;
      }
      OpNN(const OpNN &other) : OpBase<SCALAR>(other) {}

    // assignment operator
    public:
      OpNN &operator=(const OpNN &other) {
        if (this != &other)
          OpBase<SCALAR>::operator=(other);
        return *this;
      }

    // Hartree-Fock integral contributions
    void get_rhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const {
      const int s0 = this->_sites.at(0);
      const int s1 = this->_sites.at(1);
      if (s0 % 2 == 0 && s1 % 2 == 0) {
        throw pcps::Exception("both alpha case for opNN::get_rhf_ints is not ready yet, please verify code is correct");
        const int t0 = s0 / 2;
        const int t1 = s1 / 2;
        two_int[ n*(n*(n*t1+t1)+t0)+t0 ] += this->_coeff;
        two_int[ n*(n*(n*t1+t0)+t1)+t0 ] -= this->_coeff;
      } else if (s0 % 2 == 0 && s1 % 2 == 1) {
        const int t0 = s0 / 2;
        const int t1 = (s1 - 1) / 2;
        two_int[ n*(n*(n*t1+t1)+t0)+t0 ] += 0.5 * this->_coeff;
        two_int[ n*(n*(n*t0+t0)+t1)+t1 ] += 0.5 * this->_coeff;
      } else if (s0 % 2 == 1 && s1 % 2 == 0) {
        const int t0 = (s0 - 1) / 2;
        const int t1 = s1 / 2;
        two_int[ n*(n*(n*t1+t1)+t0)+t0 ] += 0.5 * this->_coeff;
        two_int[ n*(n*(n*t0+t0)+t1)+t1 ] += 0.5 * this->_coeff;
      }
    }
    void get_uhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const {
      const int s0 = this->_sites.at(0);
      const int s1 = this->_sites.at(1);
      two_int[ n*(n*(n*s1+s1)+s0)+s0 ] += 0.5 * this->_coeff;
      two_int[ n*(n*(n*s0+s0)+s1)+s1 ] += 0.5 * this->_coeff;
      two_int[ n*(n*(n*s1+s0)+s1)+s0 ] -= 0.5 * this->_coeff;
      two_int[ n*(n*(n*s0+s1)+s0)+s1 ] -= 0.5 * this->_coeff;
    }

    // function to compute the value when the operator is acted on a ket
    double action_value(const bool particles_are_spins,
                        const bool particles_are_fermions,
                        const int configs_per_site,
                        const int * occs) const {
      return occs[0] * occs[1];
      //return (occs[0] == occs[1] ? 0.25 : -0.25);
    }

  };

  //--------------------------------------------------------------------------------
  // pcps::OpRL -- class for a pair of raising and lowering operators.
  //               The first index corresponds to the raising operator and the
  //               second index corresponds to the lowering operator.
  //               For example, pcps::OpRL(2,5,1.0) would create the operator
  //               1.0 C^(+)_2 C_5
  //--------------------------------------------------------------------------------

  template<class SCALAR> class OpRL : public OpBase<SCALAR> {

    // constructors
    public:
      OpRL() : OpBase<SCALAR>(2, 0.0) {}
      OpRL(const int site0, const int site1, const SCALAR coeff)
        : OpBase<SCALAR>(2, coeff)
      {
        if (site0 == site1)
          throw pcps::Exception("The two sites of OpRL must be distinct");
        this->_sites.at(0) = site0;
        this->_sites.at(1) = site1;
        this->_offsets.at(0) =  1;
        this->_offsets.at(1) = -1;
      }
      OpRL(const OpRL &other) : OpBase<SCALAR>(other) {}

    // assignment operator
    public:
      OpRL &operator=(const OpRL &other) {
        if (this != &other)
          OpBase<SCALAR>::operator=(other);
        return *this;
      }

    // Hartree-Fock integral contributions
    void get_rhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const {
      const int s0 = this->_sites.at(0);
      const int s1 = this->_sites.at(1);
      if (s0 % 2 == 0 && s1 % 2 == 0) {
        one_int[ (s0 * n + s1) / 2 ] += this->_coeff;
      } else if (s0 % 2 != 1 || s1 % 2 != 1) {
        throw pcps::Exception("did not expect any operators with both alpha and beta sites in OpRL::get_rhf_ints");
      }
    }
    void get_uhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const {
      const int s0 = this->_sites.at(0);
      const int s1 = this->_sites.at(1);
      one_int[s0 * n + s1] += this->_coeff;
    }

    // function to compute the value when the operator is acted on a ket
    double action_value(const bool particles_are_spins,
                        const bool particles_are_fermions,
                        const int configs_per_site,
                        const int * occs) const {

      if (particles_are_fermions)
        return  (1.0 - occs[0]) * occs[1];

      if (particles_are_spins) {
        const double j = (configs_per_site-1) / 2.0;
        const double m0 = occs[0] - j;
        const double m1 = occs[1] - j;
        assert(occs[0] >= 0 && occs[0] < configs_per_site);
        assert(occs[1] >= 0 && occs[1] < configs_per_site);
        assert( (j - m0) * (j + m0 + 1) * (j + m1) * (j - m1 + 1) >= 0.0 );
        return std::sqrt( (j - m0) * (j + m0 + 1) * (j + m1) * (j - m1 + 1) );
      }

      throw pcps::Exception("OpRL could not determine the particle type");
      return 0.0;

    }

  };

  //--------------------------------------------------------------------------------
  // pcps::OpRLRL -- Class for two pairs of raising and lowering operators.
  //                 The 1st index corresponds to the 1st raising operator.
  //                 The 2nd index corresponds to the 1st lowering operator.
  //                 The 3rd index corresponds to the 2nd raising operator.
  //                 The 4th index corresponds to the 2nd lowering operator.
  //                 For example, pcps::OpRLRL(2,5,4,6,1.0) would create the operator
  //                 1.0 C^(+)_2 C_5 C^(+)_4 C_6
  //--------------------------------------------------------------------------------

  template<class SCALAR> class OpRLRL : public OpBase<SCALAR> {

    // constructors
    public:
      OpRLRL() : OpBase<SCALAR>(4, 0.0) {}
      OpRLRL(const int site0, const int site1, const int site2, const int site3, const SCALAR coeff)
        : OpBase<SCALAR>(4, coeff)
      {
        this->_sites.at(0) = site0;
        this->_sites.at(1) = site1;
        this->_sites.at(2) = site2;
        this->_sites.at(3) = site3;
        for (int i = 0; i < 4; i++)
          for (int j = i+1; j < 4; j++)
            if ( this->_sites.at(i) == this->_sites.at(j) )
              throw pcps::Exception("The four sites of OpRLRL must be distinct");
        this->_offsets.at(0) =  1;
        this->_offsets.at(1) = -1;
        this->_offsets.at(2) =  1;
        this->_offsets.at(3) = -1;
      }
      OpRLRL(const OpRLRL &other) : OpBase<SCALAR>(other) {}

    // assignment operator
    public:
      OpRLRL &operator=(const OpRLRL &other) {
        if (this != &other)
          OpBase<SCALAR>::operator=(other);
        return *this;
      }

    // Hartree-Fock integral contributions
    void get_rhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const {
      throw pcps::Exception("OpRLRL::get_rhf_ints not yet implemented");
    }
    void get_uhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const {
      throw pcps::Exception("OpRLRL::get_uhf_ints not yet implemented");
    }

    // function to compute the value when the operator is acted on a ket
    double action_value(const bool particles_are_spins,
                        const bool particles_are_fermions,
                        const int configs_per_site,
                        const int * occs) const {
      throw pcps::Exception("OpRLRL::action_value not yet implemented");
      return 0.0;
    }

  };

  //--------------------------------------------------------------------------------
  // pcps::OpIdentity -- class for the identity operator
  //--------------------------------------------------------------------------------

  template<class SCALAR> class OpIdentity : public OpBase<SCALAR> {

    // constructors
    public:
      OpIdentity(const SCALAR coeff) : OpBase<SCALAR>(0, coeff) {}
      OpIdentity(const OpIdentity &other) : OpBase<SCALAR>(other) {}

    // assignment operator
    public:
      OpIdentity &operator=(const OpIdentity &other) {
        if (this != &other)
          OpBase<SCALAR>::operator=(other);
        return *this;
      }

    // Hartree-Fock integral contributions
    void get_rhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const {
      throw pcps::Exception("no implementation for OpIdentity::get_rhf_ints");
    }
    void get_uhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const {
      throw pcps::Exception("no implementation for OpIdentity::get_uhf_ints");
    }

    // function to compute the value when the operator is acted on a ket
    double action_value(const bool particles_are_spins,
                        const bool particles_are_fermions,
                        const int configs_per_site,
                        const int * occs) const {
      return 1.0;
    }

  };

  //--------------------------------------------------------------------------------
  // pcps::OpN -- class for the number operator
  //--------------------------------------------------------------------------------

  template<class SCALAR> class OpN : public OpBase<SCALAR> {

    // constructors
    public:
      OpN() : OpBase<SCALAR>(1, 0.0) {}
      OpN(const int site0, const SCALAR coeff)
        : OpBase<SCALAR>(1, coeff)
      {
        this->_sites.at(0) = site0;
        this->_offsets.at(0) = 0;
      }
      OpN(const OpN &other) : OpBase<SCALAR>(other) {}

    // assignment operator
    public:
      OpN &operator=(const OpN &other) {
        if (this != &other)
          OpBase<SCALAR>::operator=(other);
        return *this;
      }

    // Hartree-Fock integral contributions
    void get_rhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const {
      throw pcps::Exception("no implementation for OpN::get_rhf_ints");
    }
    void get_uhf_ints(const int n, std::vector<SCALAR> & one_int, std::vector<SCALAR> & two_int) const {
      throw pcps::Exception("no implementation for OpN::get_uhf_ints");
    }

    // function to compute the value when the operator is acted on a ket
    double action_value(const bool particles_are_spins,
                        const bool particles_are_fermions,
                        const int configs_per_site,
                        const int * occs) const {
      return occs[0];
    }

  };

  //-------------------------------------------------------------------------------
  // action_value_object -- Computes the action value of the given operator when
  //                        acting on a ket.  
  //                        This struct is templated on the operator type (OP)
  //                        and the particle type (PART).
  //                        If there is no specialization for a particular
  //                        combination of OP and PART, then an error is raised.
  //-------------------------------------------------------------------------------

  //-------------------------------------------------------------------------------
  // action_value_object -- raise an error if no specialization is found
  //-------------------------------------------------------------------------------

  template<class SCALAR, class OP, class PART> class action_value_object {

    public:
      action_value_object(const int configs_per_site) {}
      double get_value(const int * const op_site_ket_occs, const int configs_per_site) const {
        throw pcps::Exception( (boost::format("no implementation for action_value_object with template parameters OP = %s, PART = %s")
                                 % pcps::type_lookup(typeid(OP)) % pcps::type_lookup(typeid(PART)) ).str() );
        return 0.0;
      }

  };

  //-------------------------------------------------------------------------------
  // action_value_object -- specialization for OpIdentity
  //-------------------------------------------------------------------------------

  template<class SCALAR, class PART> class action_value_object<SCALAR, pcps::OpIdentity<SCALAR>, PART> {

    public:
      action_value_object(const int configs_per_site) {}
      double get_value(const int * const op_site_ket_occs, const int configs_per_site) const {
        return 1.0;
      }

  };

  //-------------------------------------------------------------------------------
  // action_value_object -- specialization for OpNN
  //-------------------------------------------------------------------------------

  template<class SCALAR, class PART> class action_value_object<SCALAR, pcps::OpNN<SCALAR>, PART> {

    public:
      action_value_object(const int configs_per_site) {}
      double get_value(const int * const op_site_ket_occs, const int configs_per_site) const {
        return double( op_site_ket_occs[0] * op_site_ket_occs[1] );
      }

  };

  //-------------------------------------------------------------------------------
  // action_value_object -- specialization for OpN
  //-------------------------------------------------------------------------------

  template<class SCALAR, class PART> class action_value_object<SCALAR, pcps::OpN<SCALAR>, PART> {

    public:
      action_value_object(const int configs_per_site) {}
      double get_value(const int * const op_site_ket_occs, const int configs_per_site) const {
        return double( op_site_ket_occs[0] );
      }

  };

  //-------------------------------------------------------------------------------
  // action_value_object -- specialization for OpRL and fermions
  //-------------------------------------------------------------------------------

  template<class SCALAR> class action_value_object<SCALAR, pcps::OpRL<SCALAR>, pcps::Part_Fermion> {

    public:
      action_value_object(const int configs_per_site) {}
      double get_value(const int * const op_site_ket_occs, const int configs_per_site) const {
        return (1.0 - op_site_ket_occs[0]) * op_site_ket_occs[1];
      }

  };

  //-------------------------------------------------------------------------------
  // action_value_object -- specialization for OpRL and spinless fermions
  //-------------------------------------------------------------------------------

  template<class SCALAR> class action_value_object<SCALAR, pcps::OpRL<SCALAR>, pcps::Part_SpinlessFermion> {

    public:
      action_value_object(const int configs_per_site) {}
      double get_value(const int * const op_site_ket_occs, const int configs_per_site) const {
        return (1.0 - op_site_ket_occs[0]) * op_site_ket_occs[1];
      }

  };

  //-------------------------------------------------------------------------------
  // action_value_object -- specialization for OpRL and spins
  //-------------------------------------------------------------------------------

  template<class SCALAR> class action_value_object<SCALAR, pcps::OpRL<SCALAR>, pcps::Part_Spin> {

    public:
      action_value_object(const int configs_per_site)
        : _values(configs_per_site*configs_per_site, 0.0)
      {
        const double j = (configs_per_site-1) / 2.0;
        for (int occ0 = 0; occ0 < configs_per_site; occ0++)
        for (int occ1 = 0; occ1 < configs_per_site; occ1++) {
          const double m0 = occ0 - j;
          const double m1 = occ1 - j;
          _values[ occ0 * configs_per_site + occ1 ] = std::sqrt( (j - m0) * (j + m0 + 1) * (j + m1) * (j - m1 + 1) );
        }
      }
      double get_value(const int * const op_site_ket_occs, const int configs_per_site) const {
        //const double j = (configs_per_site-1) / 2.0;
        //const double m0 = op_site_ket_occs[0] - j;
        //const double m1 = op_site_ket_occs[1] - j;
        //return std::sqrt( (j - m0) * (j + m0 + 1) * (j + m1) * (j - m1 + 1) );
        return _values[ op_site_ket_occs[0] * configs_per_site + op_site_ket_occs[1] ];
      }
    private:
      std::vector<double> _values;

  };

  //-------------------------------------------------------------------------------
  // action_value_object -- specialization for OpRLRL and fermions
  //-------------------------------------------------------------------------------

  template<class SCALAR> class action_value_object<SCALAR, pcps::OpRLRL<SCALAR>, pcps::Part_Fermion> {

    public:
      action_value_object(const int configs_per_site) {}
      double get_value(const int * const op_site_ket_occs, const int configs_per_site) const {
        return (1.0 - op_site_ket_occs[0]) * op_site_ket_occs[1] * (1.0 - op_site_ket_occs[2]) * op_site_ket_occs[3];
      }

  };

  //-------------------------------------------------------------------------------
  // pcps::send_operator -- send an operator to another process
  //-------------------------------------------------------------------------------

  template <class SCALAR> inline void send_operator(const MPI::Comm & comm,
                                                    const typename pcps::OpBase<SCALAR>::const_ptr & op,
                                                    const int dest,
                                                    const int tag) {

    // get the numeric code for the operator's type
    int type_code;
    if      ( typeid(*op) == typeid(pcps::OpNN<SCALAR>)       ) type_code = 1;
    else if ( typeid(*op) == typeid(pcps::OpRL<SCALAR>)       ) type_code = 2;
    else if ( typeid(*op) == typeid(pcps::OpIdentity<SCALAR>) ) type_code = 3;
    else if ( typeid(*op) == typeid(pcps::OpN<SCALAR>)        ) type_code = 4;
    else if ( typeid(*op) == typeid(pcps::OpRLRL<SCALAR>)     ) type_code = 5;
    else throw pcps::Exception("unknow operator type in pcps::send_operator");

    // send the type code
    comm.Send((void *)(&type_code), 1, MPI::INT, dest, tag);

    // send the operator's coefficient
    SCALAR coeff = op->coeff();
    comm.Send((void *)(&coeff), sizeof(SCALAR), MPI::CHAR, dest, tag);

    // send the operator's sites
    pcps::send_vec(comm, op->sites(), dest, tag);

  }

  //-------------------------------------------------------------------------------
  // pcps::recv_operator -- receive an operator from another process
  //-------------------------------------------------------------------------------

  template<class SCALAR> inline void recv_operator(const MPI::Comm & comm,
                                                   typename pcps::OpBase<SCALAR>::const_ptr & op,
                                                   const int source,
                                                   const int tag) {

    // get the numeric code for the operator's type
    int type_code;
    comm.Recv((void *)(&type_code), 1, MPI::INT, source, tag);

    // get the operator's coefficient
    SCALAR coeff;
    comm.Recv((void *)(&coeff), sizeof(SCALAR), MPI::CHAR, source, tag);

    // get the operator's sites
    std::vector<int> sites;
    pcps::recv_vec(comm, sites, source, tag);

    // create the operator (remember to deallocate it later)
    typedef typename pcps::OpBase<SCALAR>::const_ptr BasePtr;
    if      (type_code == 1) { BasePtr temp_ptr(new pcps::OpNN<SCALAR>(sites.at(0), sites.at(1), coeff)                          ); op = temp_ptr; }
    else if (type_code == 2) { BasePtr temp_ptr(new pcps::OpRL<SCALAR>(sites.at(0), sites.at(1), coeff)                          ); op = temp_ptr; }
    else if (type_code == 3) { BasePtr temp_ptr(new pcps::OpIdentity<SCALAR>(coeff)                                              ); op = temp_ptr; }
    else if (type_code == 4) { BasePtr temp_ptr(new pcps::OpN<SCALAR>(sites.at(0), coeff)                                        ); op = temp_ptr; }
    else if (type_code == 5) { BasePtr temp_ptr(new pcps::OpRLRL<SCALAR>(sites.at(0),sites.at(1),sites.at(2),sites.at(3), coeff) ); op = temp_ptr; }
    else throw pcps::Exception("unknow type code in pcps::recv_operator");

  }

}

#endif
